Persistence and Quiescence of Seismicity on Fault Systems 
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We study the statistics of simulated earthquakes in a quasistatic model of two parallel heterogeneous 
faults within a slowly driven elastic tectonic plate. The probability that one fault remains dormant 
while the other is active for a time At following the previous activity shift is proportional to At~'^^^', 
a result that is robust in the presence of annealed noise and strength weakening. A mean field theory 
accounts for the observed dependence of the persistence exponent a; as a function of heterogeneity 
and distance between faults. These results continue to hold if the number of competing faults is 
increased. 
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Each of the largest earthquakes in recent decades in 
Southern Cahfornia has occurred on a different fault of a 
complex network. In contrast, most efforts at modeling 
seismicity and its complexity have been directed at stud- 
ies of self-organization on a single fault. We show quan- 
titatively that coupling between inhomogeneous faults 
leads to persistent complementary alternation of quies- 
cence and activity between the faults. This is related to 
the persistence phenomenon discovered in a large vari- 
ety of systems which specifies how long a relaxing 
dynamical system remains in a neighborhood of its ini- 
tial configuration. Our persistence exponent is found to 
vary as a function of heterogeneity and distance between 
faults, thus defining a novel universality class. 

Competition between faults has been previously ob- 
served qualitatively in a quasi-static self-organizing 
earthquake- fault model and in a dynamic model of 
the seismicity on two coupled parallel faults Our 
present model is a hybrid of these models of interac- 
tive faults. As in the driving stress is applied and 
the stress upon fracture is redistributed quasistatically 
through an elastic medium of large extent; as in rup- 
tures can occur on either of two predefined coupled, par- 
allel faults, each with variable friction. Our model is a 
crack version Q of the 2-D quasistatic earthquake model 
of a tectonic plate with long range elastic forces A 
thin tectonic plate located in the x-y plane undergoes 
antiplane scalar shear deformation along the z-axis. The 
plate is discretized in L x L plaquettes of lattice size a 
that are oriented at 45 degrees with respect to the edges; 
elastic bonds connect the nodes of the lattice. We im- 
pose a slow, uniform velocity Vz between the two oppo- 
site edges parallel to the a;-axis. The system is periodic 
in the x-direction. We imbed two linear faults parallel to 
the X-axis; we can consider more than two faults. The 
faults have finite, spatially variable, rupture strengths; 
the other parts of the system are assumed to be much 



stronger and never rupture. To ensure isotropy, we rup- 
ture nodes at the intersections of bonds, rather than indi- 
vidual bonds. The system is loaded until the total stress 
on one fault node reaches a predefined threshold. At 
this instant, loading is suspended and the correspond- 
ing bonds rupture with a complete relaxation of stress 
followed by a redistribution of stress within the lattice. 
If additional nodes become critical, the fracture will con- 
tinue to grow until no further ruptures occur; this defines 
the end of an event. At this time, the broken bonds are 
restored to full strength, and a constant, externally de- 
rived, irreversible slip is applied to each of the broken 
bonds to correspond to the force drop associated with 
the event. Loading now resumes and the process con- 
tinues. All bonds have the same elastic constant in all 
simulations. 

In contrast to earlier simulations |2|jJ] , there is no spa- 
tially variable quenched disorder on the faults and we 
instead introduce annealing of the friction by randomiz- 
ing the thresholds of ruptured nodes after each rupture 
to the values BR{x) where R{x) is uniformly distributed 
in the interval [1 — r, 1 + r] , and B is the average value of 
the force rupture threshold. Physically, the randomiza- 
tion accounts for the fact that a slip associated with an 
earthquake brings different asperities in contact that were 
previously separated, thereby changing the local sliding 
threshold. 

Crack models of this type, with relatively homoge- 
neous fracture thresholds, often display runaways [§,D, 
which are finite-size events that span the entire length of 
a fault, and which reset the stress field everywhere along 
the fault to zero. Runaways arise because the stress en- 
hancement a^J L/a at the crack tips of fractures of size 
L can overcome any heterogeneity bounded by B{1 -f r) 
for a sufhciently large crack. Real earthquake fractures 
have other length scales, and thus do not have force con- 
centrations that become excessively large due to scaling 
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by the length of the fracture. We are concerned that 
our finding of the persistent complementary alternation 
of quiescence and activity with a power law distribution 
of time scales might result from the strong overprint of 
the stress correlations induced by runaways, and thus 
might be dominated by finite size effects. We have thus 
extended our model by introducing strength weakening : 
as soon as the force on a bond reaches 50% of its ini- 
tial breaking strength, the force threshold for rupture of 
a bond is assumed to decrease linearly with time with 
a decay rate that is twice the loading rate. The results 
below do not depend on the specific values of these two 
parameters as long as strength weakening is present at a 
rate larger than but comparable to the tectonic loading 
rate, and the threshold for the onset of weakening is suf- 
ficiently small. Strength weakening simulates the effect 
of stress corrosion in fault zones. Weakening causes rup- 
tures to occur earlier, and thus the energy available for 
prolonged fracture growth is reduced. In our simulations, 
we have not observed runaways in these conditions. 
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FIG. 1. Space-time history of fault activation in a system 
of 24 X 24 nodes with two faults parallel to the x-axis. We 
have simulated systems up to 96 x 96 with similar results. 



Due to the conservative nature of elastic deformation, 
at least one fault must be active to release the average 
tectonic loading In the case of two competing faults, 
a spontaneous symmetry-breaking occurs, which is an 
example of a noise-induced transition ||^ : there is a well- 
defined threshold Tc for the amplitude of the annealed 
noise below which all rupture events take place on only 
one of the two faults and there is no activity on the other 
fault. The two faults are not decoupled, otherwise they 
would both be independently active of each other; only 
their stress fluctuations are decoupled. Tc increases with 
the distance between the two faults. Fig. 1 is an example 
of the space-time representation of the activity on two 
competing faults for r = 0.45 > Tc- Above rc, the sym- 
metry is restored and the rupture events occur on average 
equally between the two faults. 

Shortly after the time that a quiescent fault becomes 
active, activity on the other may cease along its entire 
length. The onset of quiescence does not take place in- 
stantly, but does so in a finite time after the comple- 
mentary flip to activity on the other fault. Activity 



on the quiescent fault is initiated locally at first, and 
then fractures may appear progressively along its entire 
length. Occasionally, ruptures on a quiescent fault start 
to spread, but the effort does not lead to full-fledged ac- 
tivity, as at time 455 on the first fault of the pair. Once 
an initial seed fracture has been planted on the quiescent 
fault, it can initiate rupture on its neighbors because of 
the local high stresses at the edge of the fracture. Op- 
posite these sites, fractures on the active fault are extin- 
guished because they are now in the stress shadow of the 
newly occurred fractures on the quiescent member. The 
two faults are now in competition : fractures on the for- 
merly quiescent fault attempt to spread and fractures on 
the formerly totally active fault, attempt to spread into 
the "hole" of inactivity at the sites complementary to 
the 'quiescent' fault. When the two faults are both pre- 
occupied with activity partially distributed along their 
lengths, it is impossible to decide which had been the 
active and which the quiescent fault. The formerly qui- 
escent fault may return to a state of quiescence or it may 
become completely active. Thus there is a finite interval 
of time between the triggering of activity on the quies- 
cent fault and the extinction of activity on an active fault 
in the case of a 'fiip'. In many cases, there is an outburst 
of activity, but no fiip takes place. 

100000 H ' ' ^ ' • — ^ ' • ^ ' ' ^ • ' ^: 



10000 - 
I 1000 - 




0.1 1 10 100 1000 10000 

Time Interval (B/y,k) 

FIG. 2. Cumulative distribution of quiescent time intervals 
measured in units of B /Vzk along one of the faults. The 
distribution for the other fault is similar. 

To define fault activity, we take a moving time average 
of the total slip in all rupture events on a single fault. A 
fault is defined to be active when the window-average slip 
is larger than a given threshold. Fig. 2 shows the cumula- 
tive distribution of quiescent time intervals At along one 
of the faults. The results are robust with respect to the 
threshold and size of the window, for large enough At in 
comparison with the width of the window. For small At, 
the distribution is that of the time intervals between two 
consecutive events on the same fault, which are scaled 
by the tectonic loading time B/V^k, where k is the elas- 
tic constant. Greater than this time, the distribution is 
exclusively that for intervals of flipping of activity be- 
tween the two faults. The differential distribution of the 
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quiescent intervals At is a power law 

P{At) ~ l/Ati+^('-'<*) , (1) 

for all r > rdd). The larger the value of r, the greater 
the number of flips per unit time. The exponent x thus 
increases with r (see Fig. 3), decreases with increasing 
fault separation d and is independent of the level of the 
relaxation threshold and the relaxation rate within the 
error bars. 
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FIG. 3. Dependence of the exponent x on fault heterogene- 
ity r for a fault separation of 0.062 of the fault length. The 
maximum likelihood exponent is plotted for both faults. 

Tectonic loading plays no role in the flipping of a fault 
from a dormant to an active state since accumulating tec- 
tonic stress is continuously released on the active fault. 
An active fault may flip to quiescence shortly after the 
time that the dormant fault receives a large enough force 
from the active one to trigger one site on it into activ- 
ity. We focus attention on the fluctuations in force on the 
dormant fault. For long fault lengths, the events that are 
remotely separated on a single fault are relatively inde- 
pendent. From the central-limit theorem, we expect the 
distribution of fluctuations in stress on the active fault to 
be gaussian and this is indeed the case. As the annealing 
parameter r increases, the increase in the range of rup- 
ture thresholds translates to an increase in the variance 
of the distribution of stresses; the larger the fluctuations 
in stress on the active fault, the greater the chance of ini- 
tiating rupture on the inactive fault. For small enough 
r, the activity on a single fault will runaway all too eas- 
ily. For larger but also low values of r, the variance of 
the stress on the active fault is too small to have signifi- 
cant fluctuations of the stress field to trigger activity on 
the quiescent fault. This leads to the following predic- 
tion: for two faults with the same average threshold B 
but with different r's, we expect the more heterogeneous 
member, that is the one with the largest r, to produce 
greater stress fiuctuations on the active fault that might 
activate the quiescent fault. Therefore, the fault with the 
greater heterogeneity must be the least active of the two. 
As a test of this assertion, we find that only fault 2 is 
active in the case ri = 75% and r2 = 20%. Both faults 



are active for n = 75% and r2 = 55%, but fault 2 is 
active for more than 90% of the time. As expected, the 
exponent x decreases with increasing distance between 
the faults. 

The amplitudes of these fluctuations in stress depend 

on both the length and position of the events on the 
active fault. Reactivation of the quiescent fault will occur 
at the rare time that a large fracture occurs on the active 
fault (but not so large as to cause runaway) so that a 
large force will be generated near the edge of the crack; 
the edge of this fracture has about the same x-coordinate 
as the site with low quenched strength on the quiescent 
fault. Thus, the quiescent fault will become active when 
a stress fluctuation aa induced by the active fault on it 
becomes larger than (Bq (x) — Uq (x, t))„iin , where o"g {x, t) 
is the stress quenched on the quiescent fault. We use a 
mean field theory to describe the state of stress (Ja{x,t) 
due to an active fault on the opposing, dormant fault. 
Suppose CTq fluctuates randomly due to fault activity with 

the distribution D„ = C ' ' for < aa < 

y/2'KSa 

G{d)B{l+r). C is the normalizing factor which accounts 

for the compactness of the distribution support and G{d) 
is the Green function that transfers stress from the active 
fault to the quiescent one. ctq is the mean stress equal 
to BG{d). The standard deviation must perforce be 
much less than Br. There is a quenched stress 'barrier' 
E = Bq{x) — (Tq{x) ou the quiescent fault, where Bq and 
tJq are the frozen fracture thresholds and stresses on it, 
the latter due to the nonuniform slip on it at the time of 
onset of quiescence. We model the distribution of barriers 
S as £>s(S) = (g/ e-(^-^)'/^'^' ^ijgj.g ^ Tijg 

variances and s|; are in general different because they 
reflect different stress fluctuations: is averaged over 
all possible stress conflgurations induced by the active 
fault on the dormant fault while s|. is controlled by the 
variability of the quenched stress left over from the last 
events before quiescence. 

The distribution of quiescent intervals At, for some 
fixed S, is the same as the distribution of the number 
(proportional to At) of events on the active fault that 
are necessary to get a larger than some fixed S: 

P^{At)=l dE' D^{E'){j D^{E')dE'\ 

(2) 

The complete distribution of times to reactivation by an 
isolated event on the quiescent fault is the sum of (2) 
over all possible values of S weighted by its corresponding 
distribution: 

P{At) = [ £'s(S) P^{At) . (3) 

Jo 

For At large, (2) and (3) give (1) with 
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For the WeibuU distribution = e v y ^ and a sim- 
ilar expression for we get a; = s"/ss". Very long 
tailed power law distributions for P(S) gives an expo- 
nential distribution for (1). 

Since Sa is an increasing function of r, we expect from 
(4) that X will also increase with r for smaller values of 
r, as observed (Fig. 3). For r > 0.4, the fluctuations in 
stress saturate and x becomes independent of r. Also Sa 
is a decreasing function of fault separation (d) due to the 
smoothing of stress inhomogeneities with distance. Our 
results are robust whether the breaking strength distri- 
bution is power law or uniform, and the distributions 
need not be the same from fault to fault. Our results are 
found to hold for up to 10 faults in competition which is 
the maximum number of faults we have considered. 

The number of flips m{t) during a time interval t in- 
creases as m{t) ^ t", and hence is a self-affine function of 
dimension x < I on the time axis. To see this, consider m 
successive events separated in time by Ati,i — 1, ...,m, 
where Ati + At2 + ... + Ai,„ — t — m{At), and 



{At) 



dt 



t 



Atl 



among m trials is typically 



t^, which is valid for x < 1. 



Since the maximum Atmax 

1 TOO rif' 

given by mj^^_^ ~ 

t = jn{At) ~ m~, i.e. m 
For X > 1, TO ~ 

The power law distribution of time intervals between 
flips implies nonstationarity and aging. Because of the 
self-similarity embodied in the power-law distributions, 
we can state that the longer since the last flip of activity, 
the longer the expected time till the next In other 
words, any expectation of a flip that is estimated today 
depends on the past in a manner which does not decay. 
This is a hallmark of aging. The mechanism is similar to 
the "weak breaking of ergodicity" in spin glasses that oc- 
curs when the exponent x of the distribution of trapping 
times in metastable states is less than one [||. 

Observations of activity flipping between fault 
branches on a time scale larger than inter-earthquake 
times are difficult to document due to the scarcity of 
reliable historical data. In a well-documented case 
localized flips in activity have been reported. Temporal 
clustering of earthquakes on some faults has been identi- 
fied (see ]lO[ | for a summary), but intermittency and fault 
interactions in the spirit of Fig. 1 have not been observed, 
which would be needed to establish the relevance of our 
flipping mechanism. 

We have shown that the competition between faults 
leads to the persistence phenomenon with a power law 
distribution of the durations of quiescent phases. Our re- 
sults suggest that complexity in earthquakes and faulting 
may be a generic outcome of the dynamics of interacting 
faults. We believe that a broad study of paleoseismicity 



on many faults is needed to identify intermittency as we 
have described it. 
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